#ifndef __MGGS_Method__
#define __MGGS_Method__

#include "EquationSolver.h"
#include <cmath>
#define PI acos(-1)

template<typename T>
class Bisection: public EquationSolver<T>{
private:
  T left;
  T right;
  T delta;
public:
  Bisection();
  Bisection(const int _max_iter, const T _tol, T (*_pf)(T), const T _left, const T _right, const T _delta);
  Bisection(const int _max_iter, const T _tol, T (*_pf)(T), T (*_pdf)(T), const T _left, const T _right, const T _delta);
  ~Bisection();
  void set_left(const T _left);
  void set_right(const T _right);
  void set_delta(const T _delta);
  T get_left() const;
  T get_right() const;
  T get_delta() const;
  bool Is_proper() const;
  T solve(); 
};


template<typename T>
Bisection<T>::Bisection():EquationSolver<T>(){
  left = 0;
  right = 0;
  delta = 0;
}

template<typename T>
Bisection<T>::Bisection(const int _max_iter, const T _tol, T (*_pf)(T), const T _left, const T
                        _right, const T _delta):EquationSolver<T>(_max_iter, _tol, _pf){
  left = _left;
  right = _right;
  delta = _delta;
}

template<typename T>
Bisection<T>::Bisection(const int _max_iter, const T _tol, T (*_pf)(T), T (*_pdf)(T), const T _left, const T _right, const T _delta):EquationSolver<T>( _max_iter, _tol, _pf, _pdf){
  left = _left;
  right = _right;
  delta = _delta;
}

template<typename T>
Bisection<T>::~Bisection(){};

template<typename T>
void Bisection<T>::set_left(const T _left){
  left = _left;
}

template<typename T>
void Bisection<T>::set_right(const T _right){
  right = _right;
}

template<typename T>
void Bisection<T>::set_delta(const T _delta){
  delta = _delta;
}


template<typename T>
T Bisection<T>::get_left() const{
  return left;
}

template<typename T>
T Bisection<T>::get_right() const{
  return right;
}

template<typename T>
T Bisection<T>::get_delta() const{
  return delta;
}

template<typename T>
bool Bisection<T>::Is_proper() const{
  if (left >= right)
    return false;
  T (*pf)(T) = this->get_function();
  T left_value = (*pf)(left);
  T right_value = (*pf)(right);
  if ((left_value > 0 && right_value < 0 ) || (left_value < 0 && right_value > 0))
    return true;
  else
    return false;
}

template<typename T>
T Bisection<T>::solve(){
    if (this->Is_proper() == false){
    cout << "Error! The question is not proper." << endl;
    return 0;
  }
  T (*pf)(T) = this->get_function();
  T left_value = (*pf)(left);
  T right_value = (*pf)(right);
  T mean = 0;
  T mean_value = 0;
  T s_left = left;
  T s_right = right;
  for (int i = 1; i <= EquationSolver<T>::max_iter ; i++){
    T h = right - left;
    mean = left + 0.5*h;
    mean_value = (*pf)(mean);
    if (h < delta || fabs(mean_value) < EquationSolver<T>::tol)
      break;
    else if ((left_value > 0 && mean_value < 0 ) || (left_value < 0 && mean_value > 0)){
      right = mean;
      right_value = mean_value;
    }
    else{
      left = mean;
      left_value = mean_value;
    }
  }
  left = s_left;
  right = s_right;
  if (fabs(mean_value) < EquationSolver<T>::tol)
    return mean;
  else {
    cout << "Warning! The specified accuracy is not reached." << endl;
    return mean;
  }
}


template<typename T>
class Newton: public EquationSolver<T>{
private:
  T init;
public:
  Newton();
  Newton(const int _max_iter, const T _tol, T (*_pf)(T), T (*_pdf)(T), const T _init);
  ~Newton();
  void set_init(const T _init);
  T get_init() const;
  T solve(); 
};


template<typename T>
Newton<T>::Newton():EquationSolver<T>(){
  init = 0;
}


template<typename T>
Newton<T>::Newton(const int _max_iter, const T _tol, T (*_pf)(T), T (*_pdf)(T), const T _init):EquationSolver<T>( _max_iter, _tol, _pf, _pdf){
  init = _init;
}

template<typename T>
Newton<T>::~Newton(){};


template<typename T>
void Newton<T>::set_init(const T _init){
  init = _init;
}


template<typename T>
T Newton<T>::get_init() const{
  return init;
}


template<typename T>
T Newton<T>::solve(){
  T (*pf)(T) = this->get_function();
  T (*pdf)(T) = this->get_dfunction();
  T x = init;
  T value = 0;
  T s_init = init;
  for (int i = 1; i <= EquationSolver<T>::max_iter ; i++){
    value = (*pf)(x);
    if (fabs(value) < EquationSolver<T>::tol)
      break;
    x = x - value/(*pdf)(x);
  }
  init = s_init;
  if (fabs(value) < EquationSolver<T>::tol)
    return x;
  else {
    cout << "Warning! The specified accuracy is not reached." << endl;
    return x;
  }
}


template<typename T>
class Secant: public EquationSolver<T>{
private:
  T init0;
  T init1;
  T delta;
public:
  Secant();
  Secant(const int _max_iter, const T _tol, T (*_pf)(T), const T _init0, const T _init1, const T _delta);
  Secant(const int _max_iter, const T _tol, T (*_pf)(T), T (*_pdf)(T), const T _init0, const T _init1, const T _delta);
  ~Secant();
  void set_init0(const T _init0);
  void set_init1(const T _init1);
  void set_delta(const T _delta);
  T get_init0() const;
  T get_init1() const;
  T get_delta() const;
  T solve(); 
};


template<typename T>
Secant<T>::Secant():EquationSolver<T>(){
  init0 = 0;
  init1 = 0;
  delta = 0;
}

template<typename T>
Secant<T>::Secant(const int _max_iter, const T _tol, T (*_pf)(T), const T _init0, const T
                        _init1, const T _delta):EquationSolver<T>(_max_iter, _tol, _pf){
  init0 = _init0;
  init1 = _init1;
  delta = _delta;
}

template<typename T>
Secant<T>::Secant(const int _max_iter, const T _tol, T (*_pf)(T), T (*_pdf)(T), const T _init0, const T _init1, const T _delta):EquationSolver<T>( _max_iter, _tol, _pf, _pdf){
  init0 = _init0;
  init1 = _init1;
  delta = _delta;
}

template<typename T>
Secant<T>::~Secant(){};

template<typename T>
void Secant<T>::set_init0(const T _init0){
  init0 = _init0;
}

template<typename T>
void Secant<T>::set_init1(const T _init1){
  init1 = _init1;
}

template<typename T>
void Secant<T>::set_delta(const T _delta){
  delta = _delta;
}


template<typename T>
T Secant<T>::get_init0() const{
  return init0;
}

template<typename T>
T Secant<T>::get_init1() const{
  return init1;
}

template<typename T>
T Secant<T>::get_delta() const{
  return delta;
}


template<typename T>
T Secant<T>::solve(){
  T (*pf)(T) = this->get_function();
  T x0 = init0;
  T x1 = init1;
  T x0_value = (*pf)(init0);
  T x1_value = (*pf)(init1);
  T s_init0 = init0;
  T s_init1 = init1;
  for (int i = 2; i <= EquationSolver<T>::max_iter ; i++){
    if (fabs(x0_value) < fabs(x1_value)){
      T tmp1 = x0;
      x0 = x1;
      x1 = tmp1;
      T tmp2 = x0_value;
      x0_value = x1_value;
      x1_value = tmp2;
    }
    T s = (x1 - x0)/(x1_value - x0_value);
    x0 = x1;
    x0_value = x1_value;
    x1 = x1 - x1_value * s;
    x1_value = (*pf)(x1);
    if (fabs(x0 - x1) < delta || fabs(x1_value) < EquationSolver<T>::tol)
      break;
  }
  init0 = s_init0;
  init1 = s_init1;
  if (fabs(x1_value) < EquationSolver<T>::tol)
    return x1;
  else {
    cout << "Warning! The specified accuracy is not reached." << endl;
    return x1;
  }
}




#else
// DO NOTHING.
#endif
